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This paper identifies a new pathology that can be found for numerical simulations of nonlinear 
conservation law systems. Many of the difficulties already identified in the literature (rarefaction shocks, 
carbuncle phenomena, slowly moving shocks, wall heating, etc) can be traced to insufficient numerical 
dissipation, and the current case is no different. However, the details of the case we study here are 
somewhat unique in that the solution which is found by the numerics is very weak and can fail to have 

1 ^ a derivative anywhere in the post-shock region. 
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^ 1 Introduction 

I CM„*...«C3 „„.„H»1 .,.e™ of PDE, .ee„ 

difficult to ascertain even in ID. With the notable exception of the random choice method of Glimm ^ 
I and its extensions, rigorous error bounds have remained largely elusive. This is not a situation where there 

^ is simply a hole in the theory but convergence troubles are never found in practice. There are in fact a 

number of well-known examples where numerical methods are known to behave poorly. Examples in a single 
space dimension include rarefaction shocks at sonic points the so-called wall heating phenomenon 3, 4^, 
and sub-linear convergence for linear waves [5j [6] . Often these difficulties are associated with discontinuous 
solutions or a lack of sufficient dissipation in the method because the nonlinear artificial dissipation inherent 
to the schemes goes to zero at certain points in the flow. We investigate here a new pathology where the 
dissipation is insufficient over a large portion of the domain. The result is convergence to a very weak 
solution which is nowhere differentiable. The exact nature of this weak solution is seen to depend heavily 
on the choice of time step. Admittedly this poor behavior could be eliminated with a simple linear artificial 
. viscosity, but the intent here is simply to indicate the kind of behavior that can be found. 

X 

^ 2 Governing equations and model problem 

Consider the one-dimensional Euler equations with ideal equation of state 

where u — [p, pu, pE,]'^ and f(u) = [pu.pv^ + p,u{pE + p)]'^ . Here p is the density, u the velocity, E the 
total energy per unit mass, and p the pressure. The total energy for the fluid is given by i? e + ^u^ where 
the equation of state is given by e = p(j-i) with 7 the ratio of specific heats. 

We investigate an impact problem with 7 = 1.4 on the domain x G [—.5, .5]. The initial conditions in 
primitive variables are 



[P,u,p] 



1.0, 2.0, i 



1.0, -2.0, i 



for a; < 
for X > 0. 



7_ 

Inflow conditions are applied at domain boundaries, and we integrate to time tf = 0.5. 
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Figure 1: Density (top left), pressure (top right), and a zoom of the density (bottom). In all plots the 
black line is the entropy satisfying solution and the red 'x' marks and corresponding line are the first-order 
Godunov approximation. 



3 Numerical results 

We approximate the solution using the first-order Godunov method with Roe's approximate Riemann 
solver [SI [9] . Note that the results do not change in any significant way if one instead uses an exact Riemann 
solver, nor if one moves to a second-order or high-resolution scheme. Discretization is performed on the 
computational mesh Xi — —0.5 + (« — 1) Ax for i = 1 . . . m and Ax — 1/(to — 1). Initial conditions are applied 
with exact states to the left and right of the origin. For the cases considered here, Xi for any i and the 
initial condition is applied as an exact conservative average of the left and right conservative states. 

Approximate solutions to the impact problem are shown in Figure [l] Shown are the density, pressure, 
and a restricted view of the density for a CFL number of 0.9 and m = 401. The eventual limiting behavior 
can already seen in the zoom of the density where the solution oscillates around the entropy satisfying 
solution. Because the artificial viscosity in the Godunov method is dependent on the velocity and because 
the exact solution has no post-shock velocity for this problem, the magnitude of the oscillations does not 
decrease if the post-shock velocity converges to zero fast enough in some sense. This appears to be the case 
and the frequency of oscillation increases in an unbounded manner as the grid resolution increases. Such an 
approximation will be correct in some average sense, but will not converge in an L2 or even Li sense. The 
solution to which the numerical approximation is converging appears to be a measure valued solution, whose 
value at a point can be drawn from a statistical distribution. Intuitively one can think that as Ax — 0, the 
approximate solution in the post-shock region is converging to a solution that lives in a band containing the 
exact solution. The size of that band is dependent on the details of the discretization, and most prominently 
the time step. Furthermore the width can be zero for certain circumstances, as shown below. 

We perform a convergence study using the discrete Li norm to judge convergence. The results are shown 
in the table of Figure [2j We can see that the density does not converge below approximately 5 x 10"'^. To 
understand the details, consider a single cell as it transitions from before to after the shock. Because the 
method is conservative and seems to be converging to some weak solution, the approximate shock location 
should be correct in the limit. However, each point travels through the shock in a slightly different manner 
which leads to the oscillations in density. To see this more clearly we can modify the time step such that 
the shock travels through a computational cell in an integral number of time steps. It is determined that 
the shock speed is w .76205. For the results in Figure [l] (with CFL = 0.9 and m = 401), the time step is 
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Figure 2: L\ errors and convergence rates for density, velocity, and pressure for CFL= 0.9. 





Figure 3: Zoom of the density where the shock travels through one computational cell in exactly 5 time steps 
(left) and exactly 5.5 steps (right). 



found to be At 7.496 x 10"''. As a result, the shock travels through each cell in approximately 4.3765 time 
steps. For the case of m = 401, a time step of Ai sa 6.5613 x 10~^ corresponding to CFL w .787, will have 
the shock traveling through each computational cell in exactly 5 time steps. Figure [3] demonstrates that 
this choice does indeed remove the post-shock density oscillations. Also in this figure we show the results 
for At = 5.965 x lO^'' corresponding to CFL « .716 where the shock travels through one cell in exactly 
5.5 time steps. Here we see that every other cell travels through the shock in the same way with the result 
that the density oscillates cell to cell between two values. The table in Figure [4] shows the results of a grid 
convergence study with CFL « .787 where the shock travels through each cell in exactly 5 time steps. Here 
the density is seen to converge well even to the finest resolution. 



4 Conclusions 

The pathology of convergence of numerical approximations for impact problems studied here is an interesting 
one. We showed that depending on the time step, the weak solution found by the Godunov method can 
be convergent only in some average sense. The density in the post-shock region is found to depend on the 
number of time steps required for a shock to cross a single computational cell. By choosing the time step 
so that the shock crosses a cell in an integral number of time steps, pointwise convergence can be obtained 
everywhere except a set with zero measure. Other values of the time step result in oscillatory behavior which 
causes the derivative to not exist everywhere behind the shock. 
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Figure 4: L\ errors and convergence rates for density, velocity, and pressure for time step chosen so that the 
shock travels through each cell in exactly 5 time steps. This corresponds to CFL .787. 
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